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Abstract 

In this commentary I utilize the general methods of the mathematical theory of heterogeneous 
populations in order to point out an omission in the analysis of the mathematical model in [Dwyer 
et al (2000), Am Nat, 156(2): 105—120], which led to the conclusion in [Elderd et al (2008) Am Nat, 
172(6):829-842] that the original model must be replaced with an alternative one because of the 
new data. I show that more thorough understanding of the underlying model allows twitching the 
model parameters to account for the observed data. 
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1 Introductory remarks 

In a series of influential papers [3 El[lli5] Greg Dwyer and his co-authors proposed a general math¬ 
ematical model to account for host heterogeneity among individual gypsy moth larva with respect 
to susceptibility to the virus. These papers combined accurate mathematical modeling, laboratory 
dose-response experiments, field transmission experiments, and observations of naturally occurring 
populations. In the initial works it was concluded that “a model incorporating host heterogeneity in 
susceptibility to the virus gives a much better fit to data on virus dynamics [...] than does a classical 
model” [4] and “our experimental estimates of virus transmission rates and levels of heterogeneity in 
susceptibility in gypsy moth populations give model dynamics that closely approximate the dynamics 
of real gypsy moth populations” [3j . However, in a more recent work [5] the original model was replaced 
with an alternative one, because it was found that “Our data show that heterogeneity in infection risk 
in this insect is so high that it leads to a stable equilibrium in the models, which is inconsistent with 
the outbreaks seen in North American gypsy moth populations” [5]. In more technical terms, it was 
observed that highly heterogeneous populations (for which the coefficient of variation c is bigger than 
one) exhibit oscillations with large amplitudes, whereas the mathematical model can produce such 
behaviors only for c < 1. 

My first goal in this short note is to show that the conclusion to refute the initial mathematical 
model is based on a small omission in the original analysis. To wit, while the mathematical details 
presented in [30] are mostly correct, they are based on the assumption that the host susceptibility 
is well approximated by the gamma distribution. Abstract methods of the theory of heterogeneous 
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populations, as presented in El E m\ with applications to the disease dynamics or in El 0 with 
applications to general population and evolutionary models, yield a precise characterization of the 
possible dynamical regimes depending on the chosen initial probability distribution. In particular, 
they show that an assumption on the initial distribution of susceptibility to follow the gamma distri¬ 
bution severely restricts possible time-dependent model solutions. To support my claim I present a 
mathematical model of a highly heterogeneous population (with c > 1) and stable oscillations. My 
second goal is to attract attention of the biological community to transparent and at the same time 
general mathematical models of heterogeneous populations, which contain the initial model by Dwyer 
et al as a particular case and whose underlying theory is very well worked out. Finally, the obser¬ 
vations I discuss here should prompt to reevaluate the data in [5] or collect more data to select the 
mathematical model that better approximates available observation. 


2 Mathematical model 


In my exposition I will mostly follow the notation from [3] with several minor changes. Let s(t, u) 
be the density of susceptible hosts having the susceptibility that is characterized by the parameter 
value v. Therefore, s(0, u) defines the initial distribution of susceptibility before the disease starts, 
and 5(f) = Jq s(t, v) dv is the total density of the host population at time t. Let P(t) be the density 
of infectious cadavers at t , r be the time between infection and death, and p be the breakdown rate 
of the cadavers on the foliage. Then the mathematical model takes the form 


_;(£,*,) = -us{t,u)P(t), 
d P f 

— (t)= / vP(t - r)s(t - T,v)dis - pP(t), 

d t J n 


( 2 . 1 ) 


where Q is the set of admissible values of v, e.g., D = [0, oo). The initial conditions are 


s(0,zz) = s 0 (u) = SqPq(v) = S 0 p( 0,i/), 

P( 0) = P 0 , 


( 2 . 2 ) 


where So = 5(0) is the total initial density of the host population, and po{v) is the initial distribution 
of the parameter v in the host population (such that j n po{v)dv = 1 and po{v) > 0 when G fi). I 
will use the notation p(t, v) for the current susceptibility distribution, which is given by 


n(f jA = = s & v ) 

S(t ) f n s(t,v) du' 

Inasmuch as of the most interest is the final epidemic size x (i.e., the proportion of the population 
that gets infected during the epidemics), it is possible to allow the time to go to infinity to obtain a 
transcendental equation for x in terms of the initial conditions and model parameters. Instead of this 
approach, that was used by Dwyer et al, I will show some intermediate steps. 

First, integration of the first equation in (12.11) with respect to v implies 


— (t) = -u(t)S(t)P(t), 

■^j j(t) = H* ~ t ) 5 '(* - T ) p ( t ~t)~ nP(t), 


(2.3) 
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where 


v(t ) = / vp(t, v) dv, 

Jn 

is the current mean of the distribution of v in the host population at the moment t. It is not constant, 
but a function of time (intuitively, it must decrease, since the infection washes out first those who 
have initially higher values of v). 

In [8] [TO] it was shown that function V(t) can be actually found if the moment generation function 
M(A) of the initial distribution is known, recall that it is defined as 

M(A) = f p 0 {v)e x "dv. (2.4) 

Jn 


It turns out that 

Ht) = lo § M ( A ) IA =q{t) , (2.5) 

where q(t) solves 

4f(i) = -m (2.6) 

with the initial condition q( 0) = 0. The full details of this derivation can be found in the references 
above. Therefore, instead of two equations in (12.11) I end up with three ODE (|2.3D . (12. 5p . ()2.6I) . which 
are equivalent to formally infinite dimensional system (El]). Furthermore, it can be shown that (12.31) . 
(12.51) . (12.61) are equivalent to 

4 )(<) = -h(s(t))p(t), 

f* (2.7) 

— (t) = h(S(t - r))P(t - t) - pP(t), 

where function h(S) is given by (M _1 is the inverse function to M) 


h(S) = S 0 -^M-\X)\ x=m/So . 

This is the main theoretical result of [8] that says, quite surprisingly, that the dynamics of a hetero¬ 
geneous population can be described by the same number of ordinary differential equations. I remark 
that no special assumption was made so far about the initial susceptibility distribution. 

Now assume that the initial distribution is a gamma distribution with parameters m and k: 


Po{v) 


m k—l—um 

m 


It has the mean and the variance 


E(X) = — , Var(x) = A • 

m m z 

Its moment generating function is 

M(A)= (l - —) \ 

y m J 
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and therefore, using the theory outlined above, I can find that for the initial gamma distribution 
system (12.111 takes the form 


>-m) 


1 /k 


d P, . k 

-T rW = — 
dr m 


S(t-T)^ 1/k 


So 


( 2 . 8 ) 


S(t - r)P(t - r) - pP{t), 


which coincides with equations (3)-(4) in [5]- To reiterate: system (3)-(4) in [5J is equivalent to 
the original system m if and only if the initial distribution of v is a gamma distribution. In 
[3] it was stated that “Our results, however, can also be derived without this assumption”, where 
by “assumption” the initial gamma distributions is meant. This is the point, where an incomplete 
mathematical analysis led to an incorrect conclusion, which was based on the approximation that the 
coefficient of variation, defined as the fraction of the standard deviation to the mean, is held constant. 

Let me explain in more details. When populations evolve with time, their characteristics also 
evolve (the mean, variance, coefficient of variation are all changing in general), but the distribution 
itself stays often (not always) the same. For example, if the initial distribution of susceptibility is 
the gamma distribution, than the density defined above, is again the density of the gamma 

distribution, but now with the parameters (see [6] for mathematical details) 

E (X(t)) = -, Var(X(t)) = 7 - . 

y m — q(t) v w/ ( m-q(t )) 2 

Note that both the mean and variance decrease with time since q(t) < 0 and decreasing. Using 
these expressions, we see that c = \/Var(X (t)) /E(X (t)) remains constant in time. Therefore, the 
assumption to fix c in [3] was equivalent to assuming that the initial distribution is the gamma 
distribution. It can be proved that the gamma distribution is the only distribution for which the 
coefficient of variation does not change with time in this sort of mathematical models. 


3 Toy example 


Let me here give an example of a population with the initial coefficient of variation c > 1, which 
exhibit stable oscillations. 

For the following I will need the equation for the final epidemic size x 


1 — x = Ml — 


S 0 x + P 0 


which can be obtained from 

I consider a very general family of distributions, defined by its moment generating function 

m\ -| 

M( A) = exp 


~P U- 


v — A 


with parameters v > 0, m > —1 ,mp > 0. This is the co-called variance function distributions [1]. I 
find that 


E (X(t)) = p 


m 


v-q(t)J v-q{t )’ 


VarpT(t)) = E (X(t)) 


771+1 
v - q(t) 
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Therefore, to guarantee that the coefficient of variation decreases with time, I should take m < 0, p < 0. 
For example, parameters p = —3, m = —0.2, = 0.6 imply that E(AT(0)) = 1, c(0) = 1.1547 > 1, and 

both the mean and the coefficient of variation will decrease with time. 

To produce population dynamics, I, following [5], consider the system 


N t +i = XN t (l - x t ), 

Z t +i = 4>N t x t + 7 Z t , 

where Nt and Zt are the initial host and pathogen densities in generation t. such that every time in the 
final epidemic size equation I take Sq = N t , Pq = Z t . xt is the final epidemic size for the year t, A, (j>, 7 
are the model specific parameters, which I take A = 5.5, 0 = 35,7 = 0 to coincide with the values used 
to produce Fig. 1 in [5]. The result of my simulations is given in Fig. 13.11 I note that despite the 




Figure 3.1: Dynamics of the insect-pathogen models with the initial coefficient of variation exceeding 
1. For the details and parameter values see the text 

fact that the initial coefficient of variation is above 1 (the population is highly heterogeneous), we still 
observe stable oscillations of large amplitude. This is in contrast with Fig. IB in |5j, where a stable 
equilibrium for the model with c > 1 is shown. Therefore, in general, the conclusion from [5] that the 
model with the constant infection risk cannot produce the sustained oscillations if c > 1 should be 
replaced with the conclusion that the model with the constant infection risk and the initial gamma 
distribution of susceptibility cannot produce the sustained oscillations if c > 1 . 

4 Concluding remarks 

First of all, I would like to emphasize that my observations do not invalidate the conclusions in [5]. 
The important point here is that the mathematical model, originally suggested in [4] and analyzed in 
[3], is still capable of producing dynamical regimes, which were thought to be impossible to observe in 
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it (see my toy example above). Moreover, the general theory of heterogeneous populations [6l [8, [10] 
provides a convenient and flexible tool to accommodate very different dynamical behaviors. 

A straightforward test to determine whether the assumption on the initial gamma distribution in 
|5] should be replaced with something else may be performed by collecting estimates of the population 
coefficient of variation at different time moments during each epidemics. If these values do not deviate 
with time significantly, then it would support the reasonings in [3], If, however, it would be possible 
to convincingly see that the coefficient of variation declines as the epidemics proceeds, then it should 
prompt to reevaluate the role of the model with the constant infection risk, as originally formulated 
by Dwyer and his co-authors. Unfortunately, the available data do not seem to allow to perform such 
test at the moment. 

Acknowledgements: I would like to thank Bret Elderd for a profitable discussion while preparing 
this commentary. 
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